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ABSTRACT 



This report studies two aspects of tracking ballistic missiles during boost phase. 
The first part compares the performance of several nonlinear filtering algorithms in 
tracking a single target: the extended Kalman filter (EKF); the unscented Kalman filter 
(UKF); the particle filter (PF); and, the particle filter with UKF update (UPF). Meas- 
urements are range, azimuth and elevation. In the absence of measurement error, all al- 
gorithms work well except for the PF, which does not converge. With measurement noise 
(standard deviations of 10 meters and 1 degree) the EFK also performs poorly, while the 
UPF is the top performer (although it is also the most computationally intensive). 

The second part compares the extended information filter (EIF) with earlier work 
on track scoring to perform sensor/data fusion in a multi-hypothesis framework. Here we 
find that the EIF handily outperformed other fusion algorithms based on track scoring 
that we tested. 
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I. INTRODUCTION 



The objective of this study is to investigate various algorithms for tracking multi- 
ple ballistic missiles within the boost phase. The ballistic missile behavior and flight pro- 
file are simulated using a proprietary add-on MATLAB module, specifically, 
IMPULSE©. Sensor modeling has been extensively discussed in previous work[16][20] 
and these same models are used here. 

In the first part, the estimation performance of the following nonlinear filters is 
compared: the extended Kalman filter (EKF), the unscented Kalman filter (UKF), particle 
filter (PF), and the particle filter with UKF proposal (UPF). Tracking ballistic missiles 
using infrared (bearing only) and radar (range and bearing) sensors requires a nonlinear 
filtering approach. The most significant nonlinearity results from the conversion of sensor 
measurements of range, azimuth and elevation into Cartesian coordinate estimates of the 
target state. This study compares the relative performances of various non-linear filtering 
methods applied to tracking multiple ballistic missile targets using radar measurements. 
The purpose of including the EKF in this study is due to the wide acceptance and fre- 
quent application of the EKF. The reason for including the UKF in this study is due to 
the rapidly growing support for the UKF over the EKF in many applications of nonlinear 
state estimation. The motivation behind including the PF is its appliance to many nonlin- 
ear estimation problems in recent research studies. 

Many recent papers have studied the ballistic missile tracking problem. Farian, et 
al. studied the problem of tracking a ballistic object in the reentry phase using the EKF, 
UKF and PF in [1] . They found that the results favored the EKF. Bruno and Pavlov also 
applied a variation of the PF to the case of a reentering ballistic object in [2]. Both of 
these studies used only a second order coordinate system involving position and velocity, 
and did not consider acceleration. Saulson and Chang used the UKF, EKF and Central 
Difference Filter (CDC) to track a ballistic missile in flight from a variety of platforms in 
[3], Althiough they looked at the boost phase, their treatment considered a missile with 
constant acceleration. Siouris etal. [4] also performed an analysis of the incoming ballis- 
tic object using an extended interval Kalman Filter. The American Association of Scien- 
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tists looked at the Airborne laser problem in the boost phase but only used the EKF for 
tracking [5]. In addition, they did not consider tracking from the aircraft and the required 
sensor accuracy. Clemens and Chang examined the use of various nonlinear estimation 
filters, EKF, UKF, PF, GSPF (Gaussian Sum Particle Filter), for tracking a ballistic mis- 
sile during boost phase from a moving airborne platform [10]. 

In the second part, we propose a hierarchical data fusion architecture and method , 
when a number of different types of sensors are deployed in the vicinity of a ballistic 
missile launch. Sensor data fusion is the process of combining outputs from several sen- 
sors with information from other sources, such as information processing blocks, data- 
bases or knowledge bases, into one representational form. It is expected to achieve im- 
proved accuracy and more specific inferences than could be achieved by the use of a sin- 
gle sensor alone. In principle, fusion of data from several sensors provide significant ad- 
vantages over single source data. In addition to the statistical advantage gained by com- 
bining same-source data, the use of multiple types of sensors may increase the accuracy 
with which a quantity can be observed and characterized. In this study, seven active 
ground-based radar sensors are used to track the ballistic missile in the boost phase. We 
consider the Information Filter which is an efficient form to fuse the information from 
multiple sensors without information loss. To evaluate the performance of the Extended 
Information Filter (EIF) [12], we compare with the track score function proposed in Pat- 
sikas’s previous work [20] in which the calculation of a track scoring function is used to 
identify the sensor with the best track file. All simulations have been implemented using 
MATLAB. 
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II. SENSOR MODELS AND MULTIPLE HYPOTHESES TRACKING 



In order to detect the threat of a ballistic missile attack, there are a multitude of 
sensors available to provide early warning information if such an event were to occur. 
Two particular types of sensors are appropriate for the boost-phase missile tracking prob- 
lem; passive sensors such as IR, and active sensors such as radar. In this chapter, we in- 
vestigate surface-based sensors. Furthermore, the mathematical relationships, which de- 
termine the performance of each sensor, will be examined. The results will enable us to 
take the IMPULSE© missile flight information and introduce appropriate error as a 
means to model the sensed position information as reported by each sensor. In this sec- 
tion, the radar model and multiple hypotheses tracking (MHT) algorithm, originally de- 
veloped in [16], are reviewed. 



A. GENERATING BALLISTIC MISSILE PROFILES 

This research makes use of the IMPULSE© simulation tool for generating ballis- 
tic missile flight profiles which will later serve as the test data for the tracking algorithm. 
This software is a product of the National Air & Space Intelligence Center (NASIC) — an 
organization recognized as the cognizant analysts’ representative for threat-missile plat- 
forms. IMPULSE© is used by engineers and researchers to conduct missile guidance 
testing, targeting studies and further enable users to study classified ballistic missile per- 
formances. In our study, this software enables us to include flight profiles more sophisti- 
cated than simple, constant-rate, parabolic motion-models. The missiles in the simulation 
display flight performances that one would expect a real-world counterpart to exhibit. 
IMPULSE© generates missile flight performances to reflect realistic missile physical dy- 
namics, e.g., non-linear velocity and acceleration profiles due to changing mass, staging 
events, and interaction with the Earth’s atmosphere and gravitational field. Figure 2 de- 
picts six missiles, generated with the IMPULSE© software, which have been fired from 
suspected North Korean missile facilities. 
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Figure 1 . Main functional area implemented in this study to include the IMPULSE@ 
threat profile, the RF sensor model, and the multiple hypothesis tracking algorithm. 




Figure 2. Multiple missile attack directed toward the American continent. 



The files listed in Table 1 contain the IMPULSE© output of each missile’s flight 
as exported by the General Write Utility (GWU). Information regarding the use of this 
tool can be found in Appendix A. These files will serve as the input data to the sensor 
models and tracking algorithm. 



Missile Profile 


File 


Ballistic Missile 1 Upper Stage data 


FltlUStage.txt 


Ballistic Missile 1 Lower Stage data 


FltlLStage.txt 


Ballistic Missile 2 Upper Stage data 


FltlUStage.txt 


Ballistic Missile 2 Lower Stage data 


FltlLStage.txt 


Ballistic Missile 3 Upper Stage data 


FltlUStage.txt 


Ballistic Missile 3 Lower Stage data 


FltlLStage.txt 


Ballistic Missile 4 Upper Stage data 


Flt4UStage.txt 


Ballistic Missile 4 Lower Stage data 


Flt4LStage.txt 


Ballistic Missile 5 Upper Stage data 


FltlUStage.txt 


Ballistic Missile 5 Lower Stage data 


FltlLStage.txt 
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Ballistic Missile 6 Upper Stage data 


Flt6UStage.txt 


Ballistic Missile 6 Lower Stage data 


Flt6LStage.txt 



Table 1. 



Table 1. Output files from IMPULSE General Write Utility. 

Table 2 details the file format for each text file listed in Table 1. The simulation 
time is given in the first column. The latitude and longitude are represented in the second 
and third columns. The forth and fifth columns give the altitude in meters and the thrust 
magnitude in Newtons, respectively. 



Time, 

s 


Geodetic 
Latitude, rad 


Geodetic 
Longitude, rad 


Altitude, 

km 


Thurst 
Magnitude, N 


0 


0.70969 


2.2372 


0.02 


2.03E+05 


1 


0.70969 


2.2372 


0.027182 


2.03E+05 


2 


0.70969 


2.2372 


0.048894 


2.03E+05 


J 


0.70969 


2.2372 


0.085382 


2.03E+05 


4 


0.70970 


2.2372 


0.13561 


2.03E+05 


5 


0.70970 


2.2372 


0.20043 


2.04E+05 


6 


0.70970 


2.2372 


0.27988 


2.04E+05 


7 


0.70971 


2.2372 


0.37395 


2.04E+05 


8 


0.70971 


2.2372 


0.48265 


2.04E+05 


9 


0.70972 


2.2372 


0.60599 


2.04E+05 


10 


0.70973 


2.2372 


0.74394 


2.05E+05 



Table 2. Column format for each Ballistic Missile file listed in Table 1 as the output 

from the General Write Utility GUI. 

The ballistic missile positions in these files are given in the WGS84 geodetic co- 
ordinate system. Furthermore, the latitude and longitude values in these files are in units 
of radians with respect to the center of the Earth. The data can be easily converted to 
North/South, Degrees/Minutes/Seconds (N/S DD/MM/SS) and East/West, De- 
grees/Minutes/Seconds (E/W DD/MM/SS) format by using the radldeg.m function in the 
MATLAB mapping toolbox. 
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B. SENSOR MODELS 



In this study, seven radars, which are placed at different locations relative to the 
launch position, use the same configuration in order to make the analysis easier to verify. 
The single pulse, signal to noise ratio (SNR) for each radar system can be expressed as: 



S / N - 



Pfj] tX 2 <j 
( 47rf kT s BR 4 



( 0 . 1 ) 



where T s =T a + T e , 7^ is the system temperature, T e is the receiver noise temperature and 
T a is the antenna noise temperature, which in the work reported here, is assumed to be 
290K as was the assumption in [16]. P t denotes the peak transmitted power, G , is the an- 
tenna gain, a is the radar cross section of the object , r is the compressed pulse width. 
X is the wavelength, B is the bandwidth of the receiver, k is Boltzmann’s constant and R 
is the range to the target. Generally, the radar losses are also considered but are assumed 
to be zero for this study. The sensor-to-target slant range, R , is obtained in the simula- 
tion by using the MATLAB elevation function. This returns the actual target-to-sensor 
range. Due to the large unambiguous range required, the pulse repetition frequency (PRF) 
is modeled as PRF = 150Hz. Figure 1 1 shows the SNR computed at radar sensor 1 (RF1) 
as it observes the upper and lower stages of the ballistic missile in our simulation. An in- 
creasing SNR is observed as the missile passes near and relatively over the sensor during 
the boost phase of flight. This is observed in the first 100 seconds of the missile’s flight. 
Staging occurs at 65 seconds where a second SNR value — the green line — is initiated. 
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15 




Figure 3. Radar Sensor 1 SNR versus time while observing BM1. 

The SNR value observed on the booster — the red plot — remains at a nominal val- 
ue of 8 dB as it never leaves the sensor’s immediate field of view. Conversely, the upper 
stage of the missile continues its flight away from the sensor. The SNR is inversely pro- 
portional to the forth power of the range; as a result, the SNR decays rather rapidly be- 
yond 130 seconds of the simulation. As the SNR changes throughout the course of the 
missile trajectory, the values will be used in the computation of the error in the reported 
angular and range measurements. 

Table 3 lists the parameters of the radar considered in the simulation. 



Parameter 


Value 


Carrier Frequency,/ 


10 GHz 


Peak Power, P, 


1.0 MW 


Antenna Gain, G 


42 dB 


Radar operating Bandwidth, 
B 


15 MHz 


Receiver Noise Tempera- 
ture, T e 


20 



Table 3. 



Table 3. RF sensor parameters. 



8 



The radar cross section (RCS) of the missile during the flight is calculated by the 
program POFacets [17]. The missile used in the simulation is a sample unclassified mod- 
el. To determine RCS of the missile, it is assumed that the missile is a circular cylinder 
with a radius of 0.60 cm and a length of 20m. 

The angle 6 between the missile’s trajectory vector and the distance vector from 
the radar to missile, which is required in order to define the RCS function, is determined 
by the MATLAB function RFobservation by using the geometry of Figure 4. From Fig- 
ure 4, it is inferred that the angel of incidence is computed as the inner product between 
the position vector between the sensor and the missile and the tangent vector to the mis- 
sile trajectory curve. The vector of the curve is defined from the difference of the succes- 
sive points of the trajectory. The RCS determination is a 3D problem. In this simulation, 
it is assumed that the change of the angle takes place only on the x-y axes without taking 
the third dimension into account. 



'arcs*. (oall stic nissile) 



ci.ae of me trajectory 



• — 

e 



LOS 



Incident angle 



, R. distarce between target and radar 

\ 

\ 

\ 

\ 

> 

& Radar 



Figure 4. Calculation of the incident angle 6 . 



Using the above geometry, the MATLAB function RFobservation calculates the 

incident angle dfor the first radar sensor. A plot of this angle data is shown in Figure 5. 

9 




Figure 5. Calculated incident angle 9 . 

Having determined the incident angle and using the POFacets code, the RCS of 
the ballistic missile is obtained as shown in Figure 6 and 7. 

Figure 6 is a polar plot of the RCS as a function of the angles of incidence. In the 
vicinity of angle 0 = 90°, which has been calculated from the previous analysis, the val- 
ues of the RCS are around 36 dBsm. 




Figure 6. Polar plot of the Ballistic Missile RCS. 
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Figure 7. Linear plot of the Ballistic Missile RCS 

Figure 7 is a plot of the ballistic missile RCS with respect to the angle of inci- 
dence. From the above two plots, it is inferred that the RCS for an angle of incidence in 
the vicinity of 90° is equal to 36 dBsm. 

C. MULTIPLE HYPOTHESES TRACKING 

In this section, we briefly describe the MHT (Multiple Hypotheses Tracking) al- 
gorithm developed in [16]. 

1. Contacts, Targets, Scans and Associations 

To explain multiple target tracking, several terms must be established to provide 
clarity. First, the definition of a contact is introduced. This is an observation that con- 
sists of a detection — by a sensor — and its corresponding measurement. Commonly, a de- 
tection occurs when the signal-to-noise ratio of a sensor meets or exceeds a predefined 
threshold. The RF signal return from the missile was such that it enabled the sensor to 
detect the object. Furthermore, the measurement corresponding to this detection consists 
of the position of the object. In limiting the sensor responses to contacts, a restriction is 
imposed such that a signal return is of interest only when the object meets a predeter- 
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mined criterion. This allows the contact to be “flagged” for further examination. For the 
remainder of this study, the term measurement will also be used interchangeably to mean 
contact. These terms, as well as the next several terms given, are illustrated in Figure 8 
for a single target tracking problem. 



next state 




Figure 8. Single target tracking (in 2-D) illustrating measurement-to-target state predic- 
tion pairing, and next-state prediction correction. 

A target is a measurement from a previous sample time, which has been flagged 
and declared an object of interest; hence, the term, a priori target, is used. This object’s 
information, or track file, may consist of position, velocity and acceleration parameters. 
Certain characteristics of a measurement are “screened” and, if they fall within predeter- 
mined values, the measurement is declared a target and stored in a database to await fur- 
ther information update. A contact’s velocity information, for instance, may be used to 
mark the measurement for further observation. 

The term scan will be used to refer to successive time periods. Allowable meas- 
urements and targets will now be further restricted to specific time scans. Finally, the 
term association is used to describe the pairing of a measurement, in a present scan, with 
that of a target in a previous time scan. It is important to note that there must be no am- 
biguity in a target-to-measurement pairing from one time period to the next time period. 

There are two essential assumptions made in a single target tracking problem. 
First, there is one and only one target present in the observation region. Secondly, all sen- 
sor observations in successive scans are generated by the same target. Given the a priori 
target, j, as in Figure 8, the problem is approached as follows: based on prior information 

about the target, a next-state prediction is generated. This is given by 
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** + l = F k x k + <y* 

where x M is target next-state vector (prediction), F k is the known state transition matrix, 
co k is the plant-noise associated with the target and k is the time index. Next, a measure- 
ment is taken in the ensuing scan. Since the measurement is likely to originate from the 
target that caused the contact in the previous scan, a direct comparison may be made be- 
tween the new observed contact and the expected next-state prediction. The difference 
between the observed measurement and the expected prediction is referred to as the inno- 
vation 

The algorithm which computes the next-state prediction is updated with the new 
observed information. Adjustments are made to improve follow-on estimations. The al- 
gorithm eventually stabilizes such that each next-state prediction approaches the actual 
measurement. Obviously, tracking is straight-forward when there is no uncertainty about 
the origin of successive measurements and when targets do not accelerate in unpredict- 
able ways. 

On the other hand, in a multiple-target environment, for any given scan, the sen- 
sor responses — measurements — may be due to any target. Thus, each possible move- 
ment of established targets must be hypothesized in the observation space. The primary 
difficulty in multiple target tracking is correctly assigning successive contacts to a corre- 
sponding established target’s next state prediction as shown in Figure 9. This process of 
making possible assignments is referred to as measurement-to-target association. 
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Figure 9. Multiple-target tracking (2-D) scenario. An illustration of the difficulty in cor- 
rectly pairing new measurements with their state prediction so subsequent predictions 
converge on the true track. 

This is an ideal place to introduce the term association hypothesis in which each 
hypothesis attempts to provide a likely explanation as to the source of each new meas- 
urement. Hence, we have the algorithm’s namesake — MHT. 

An association hypothesis maps each measurement to a possible target’s next- 
state prediction. Furthermore, if the association of contacts were always obvious, the 
problem would simplify to independent single target problems — the recursive process of 
measurement, state prediction, observation, and, finally, update of each target. However, 
due to the unpredictable nature of each target in a multiple-target space, the associations 
are ambiguous. 



2. MHT Implementation 

The following paragraphs cover the implementation of the MHT using the terms 
and concepts that have been established in Section A. Furthermore, the strategy for solv- 
ing the association problem is presented in detail. 

- Generation of Association Hypotheses 



14 



The first step in the measurement-to-target association process is to form multiple 
feasible hypotheses. The assumptions made in this study are as follows: 

• Each target causes, at most, one measurement to be generated by the sen- 
sor. Clustering as presented in [21] is an example whereby multiple con- 
tacts with common measurements can be represented by a single contact. 

• Each measurement corresponds to only one known target. The possibility 
for sensor false alarms is not addressed in this study. 

• Sensor measurements are updated every second. 

• Targets, in this case ballistic missiles, do not deploy chaff or use electronic 
countermeasures. 

The following discussion is based upon [21] and provides the framework for the 
implementation of the MHT. Let Z k = {z km , m = 1,2,..., m k j denote the set of measure- 
ments m in scan time k. In the simulation, a single measurement z k m can be further de- 
composed such that z k m = m , y k m , £ km \ where £, y, and C, denote the position of 

the measurement of interest in a predefined coordinate system. In the study of ballistic 
missiles, a local vertical coordinate system — north, east and up — centered about the sen- 
sor's location will be used. Let Z k - jz 1 , Z 2 , Z 3 , ..., Z k j represent the cumulative set of 

measurements up to scan time k. Let Q k =|of, / = 1, 2, ..., /} represent the set of all 

measurement-to-target hypotheses at scan k. The index, i , denotes the hypothesis num- 
ber. This is the association of measurement in a current scan with a priori targets estab- 
lished in preceding scans. Furthermore, the hypothesis Q k of the Ath scan is taken as the 

joint hypothesis formed from all prior hypotheses £2*'' and the association hypothesis for 

the current data scan. Using Figure 19 as an illustration, each hypothesis comprises two 
possible association pairings. 

The probability of each feasible hypothesis may now be derived. Let P k be the 
probability of hypothesis, Q k , given measurements up through time k. To further clar- 
ify, allow o, to represent hypothesis- 1 where measurements is associated to target -1 
(implied by pairing with state estimate-]) and, within the same hypothesis, measurement- 
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2 is associated to target-2 (see Figure 9). The probability of this hypothesis being the 
correct pairing is represented by P k . Alternatively, hypothesis-2, Cl k 2 , with probability 

P 7 k corresponds to the alternate hypothesis where measurement - 1 is associated with tar- 
get-2 while measurement-2 originates from target- 1. Also, let P k ~ 1 represent the prob- 
ability Q*" 1 ; the past relationships. To successfully perform multiple hypothesis tracking 
where computational time must remain minimal, a recursive relationship must be estab- 
lished between P k and P k ~' to avoid reevaluating all past data whenever a present scan 

set, Z k , is received. 

The association probability,/’*, is a conditional probability of Q* given the set 

of measurements, Z k , and the hypothesis, Q k . In [22], Nagarajan, Chidambara, and 
Sharma (NCS) give the relationship 

P k =p(n k g \z(k), Q*) = P({Q k -\y/ h }\Z k ,Q k ) 

where y h = j m- 1, 2, 3, ..., m k j and represents the event that the mth meas- 
urement corresponds (originated from) to the y'th target as per association hypothesis if / h . 

Furthermore, g denotes the global index while h denotes the hypothesis index. By using 
Bayes’ rule, we can express the above equation as 

By setting p{o. k j = C , a constant, the expression becomes 
which can be further written as 

p‘ =i/>(Qj-')f( n |n‘-',z‘) 
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where we have again applied the Bayes’ rule. Since the events i// that comprise ij/ h are 
independent, we have 



1 mi 

^ m = 1 

The right-most term is made possible by observing that the current scan association is af- 
fected only by Z k and all past data have been included through the term Q* _I . From 

[22], by defining ) = P[y mJk | H*" 1 , Z k ) , we have 




"‘i 

/ ’W")n^(™.A.nr i ) 



m~] 



C 



(4.1) 



The term /3(m,j h ,Q k g i } represents the probability that measurement m corresponds to 
target j h as per the present scan hypothesis if/ h and the past scan hypothesis Q* -1 . 

To solve for the above probability, let J be the number of known a prior targets, 
and let ,, _/' w 2 , ..., j k _ x y have predicted next-state measurement vectors 

x kJl , x k j2 , x k ]N and corresponding innovation covariance matrices 
Z* ,, Z k l , ..., E* N , respectively, as per retained from the previous scan k- 1. When 
the new measurement z k m corresponds to a confirmed target j N whose existence is im- 
plied by the prior hypothesis Q*"', then p(jn,j h ,Cl k ~ x ^ is characterized by the probabil- 
ity density function [25] 







■’-l 






(4.2) 



where z m j is the measurement-to-target innovation, is the determinant of the in- 

novation covariance, and n denotes the degrees of freedom of the measurement (in this 
study, n = 3 as the sensor’s measurements consist of range, azimuth and elevation). 

Furthermore, the innovation is given by [23] [24] 
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( 4 . 3 ) 



c k,m,j 



= Z 



k , m 



-H.x 



j k\k-l 



and represents the difference between the observed measurement, z k m , and the pre- 
dicted — or expected — measurement H t . The observation matrix, H r is the used 

to select the elements of x . kk _ x for comparison. The innovation (or residual) covari- 
ance given by 



^ k,m,j - Hr. k,j^j,k\k-l^r, k,j + 



(4.4) 



The quantities z k m and above are obtained from the output of the extended Kal- 
man (EKF) filter as outlined in Appendix B, with the update estimate covariance, P ] klk _ x 
The noise covariance, R k , and the observation matrix H r k J are further explained in Eq- 
uations B.2 and B.3 of the same appendix. The appendix further provides a comprehen- 
sive review of the EKF equations and the motion model used to predict the next state po- 
sition x. of the target . More details are in [16]. 



In summary, this chapter described the parameter of the RF sensor at sea level 
used in the simulation and the multiple hypothesis tracking algorithm based on previous 
work [16][20]. 
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III. NON-LINEAR FILTERS FOR TRACKING 



In this section we describe four recursive estimators (filters) to be used for track- 
ing a ballistic object. The EKF, UKF, PF and PF with UKF proposal are applied in this 
study. A description of each of the filters used is taken in turn in the following sections. 

A. EXTENDED KALMAN FILTER (EKF) 

The EKF linearizes nonlinear dynamic models using partial derivatives and the 
expected state vector in order to perform Kalman Filtering with a Taylor Series approxi- 
mation of the time and/or measurement updates. As a result of the functional lineariza- 
tion, the EKF provides the optimal linear, or Linear Minimum Mean Square Error 
(LMMSE), solution. It always approximates the probability density to be Gaussian and is 
therefore not well suited to handle non-Gaussian distributions such as bi-modal or heavily 
skewed cases. To do this, the EKF uses partial derivatives of the measurement matrix h 
with respect to the state x in an attempt to linearize the non-linear dynamic observation 
equations. The derivative is performed at each time step, k, resulting in the Jacobian ma- 
trix, designated H(k). Once the partial derivative is found, standard Kalman Filtering 
procedures are used to determine the predicted state and it's covariance. More details are 
in the [ 16 ]. 

B. UNSCENTED KALMAN FILTER (UKF) 

The Unscented Kalman Filter is a recursive MMSE estimator that addresses some 
of the approximation issues of the EKF. Because the EKF only uses the first order terms 
of the Taylor series expansion of the nonlinear functions, it often introduces large errors 
in the estimated statistics of the posterior distributions of the states. This is especially 
evident when the models are highly nonlinear and the Taylor series expansion error be- 
comes significant. Unlike the EKF, the UKF does not approximate the non-linear proc- 
ess and observation models, it uses the true nonlinear models and, instead, approximates 
the distribution of the state random variable. In the UKF the state distribution is still rep- 
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resented by a Gaussian random variable (GRV), but it is specified using a minimal set of 
deterministically chosen sample points. These sample points completely capture the true 
mean and covariance of the GRV, and when propagated through the true nonlinear sys- 
tem, captures the posterior mean and covariance accurately to the 2 nd order for any 
nonlinearity, with errors only introduced in the 3 rd and higher orders [8]. 

The UKF uses the unscented transformation, or the "‘scaled unscented transforma- 
tion” implemented in this study, to generate deterministic sigma points to approximate a 
random variable’s probability distribution given a true (or assumed) nonlinear function 
and the prior mean and covariance. The sigma points are chosen based on the prior mean 
and respective columns of the matrix square root of the prior covariance. The Choleski 
factorization used to solve for the upper triangular matrix square root of the expected 
state covariance requires that the covariance matrix be positive definite. Square-root im- 
plementations of the UKF perform more efficiently and without the positive definite re- 
quirement. The UKF is similar to Particle Filters in that it generates points about the 
mean estimate, but requires less computational cost due to deterministic sampling of the 
sigma points as apposed to the Particle Filter’s random “particles”. Unlike Particle Fil- 
ters, the UKF assumes a Gaussian distribution, but it is able to approximate heavy tailed 
distributions better than the EKF. For the second-order unscented filter, the estimates of 
the mean and covariance are accurate to the 2 nd order in general and to the 3 rd order for 
Gaussian priors. The second-order unscented filter generally requires at least 2n + 1 
sigma points, while the fourth-order unscented filter requires at least 2n 2 + 1 sigma points 
and is able to estimate the mean to the 4 th order in general. A reduced sigma point im- 
plementation has only the minimum requirement that n+1 sigma points be used [9]. 
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Figure 10. Schematic diagram of the Unscented Transformations. 

To calculate the estimation of the state vector x(k+l) using the UK.F one begins 
by developing a matrix ^of 2L+1 sigma vectors each with a weight IV, where L is 

the dimension of the state vector x (in this case L=9). The sigma vectors and their 
weights are calculated by 

Z 0 =x 

Zi =x-(J(L + A)P x } )' i=l,....L (l) 

Z, =x + (Ja7M) i i=L+I,...,2L 

W 0 {m) =X/(L + y i) 

W 0 m = A/(L + A) + (\-cc 2 +j3) 

fV (m) = W t c = 1/{2(L + A)} i=i,...2L 
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where ^yj(L + Z)P x j is the ith row or column of the matrix square root of ( n x +k)P x . 

The weights are normalized. Each sigma point is propagated through the non-linear func- 
tion 



•V/ =&(*;)> i=0,...,2L (2) 

The estimated mean and covariance of t are computed as the sum of the propagated 
points 

?=*br,y, (3) 

7=0 

p , = iy,(y,-y){.y,-y) T 

/= o 

Implementation of the above algorithm with the current estimate and its covari- 
ance x(k | k) and P x ( k \ k } follows as [1] 



1. Compute sigma points C, t and their weights W t , (i=0,...,20) using (7) 

with a =0.001, /?=2, using and. P x (k\k} . The matrix square 

root in (1) is found using Cholesky factorization. 

2. The sigma points are propagated using the state equation such that 

<r,(*+ii*)=e|V,(*i*>]+e(*) (4) 

3. From these new points, compute the predicted state and covariance 
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x(k+\\k) = f i w l ,; l (k+)\k) 

1=0 

P(k + \\k) = f t W i [t l (k + \\k)-x(k + \\k)]* 

i = 0 

[£(*+i| *)-*<*+ i|*)] r + 2 (*) 



4. Predict the observation sigma points using observation equation 

C i (k + l\k) = H[C i {k\k)] 



5. The observation estimate and covariance follows 



y(k + l\k) = f i W£ l (k + \\k) 

i=0 



( 7 ) 



+ 1 1 A) = X ^ [C <* + 1 1 *) - + 1 1 *)] • 

/=0 

+ 1 ) k^j — y (jt + 1 1 J + R(k 1) 

P xy {k + \\k) = f j W l [i,{k + \\k)-x{k + \\k)}> 

1=0 

[C,- (A: + 1 1 A:) jp (A: + 1 1 /r)] ? 



6. Compute the UKF gain and updated state and state covariance similar 
to the standard Kalman Filter 

K(k + \) = P^ 

+ 1 1 k + \^ = x(k + \\ k} + K(k + l)[.y(A: + 1) — y(k + 1 1 A:^J 
P{k + \\k + \) = P(k\k + \)-K(k + \)P yy K r (k + l) 
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Note that no explicit calculations of Jacobians or Hessians are necessary to implement 
this algorithm. The UKF requires computation of a matrix square root which can be im- 
plemented directly using a Cholesky factorization of order n* / 6 . However, the covari- 
ance matrices can be expressed recursively, and thus the square-root can be computed in 
2 

order n x by performing a recursive update to the Cholesky factorization. So, not only 

does the UKF outperform the EKF in accuracy and robustness, it does so at no extra 
computational cost. 



C. PARTICLE FILTER (PF) 

The previous nonlinear Filtering algorithms rely on Gaussian approximations. In 
this section, we shall present the particle filtering method. It is a set of state estimation 
methods using a weighted set of samples (particles) in a Monte Carlo methodology to ap- 
proximate the posterior distribution of the state [10][1 1]. A large set of particles are gen- 
erated and weighted to represent the posterior distribution of the state. Each particle is 
independently propagated through time in accordance with the dynamic state space 
model and the system process noise. After a new observation y(k+l) arrives, importance 
evaluation and resampling steps are used to assess the correlation of each individual par- 
ticle to the observation. New particles are generated to replace the old ones and the parti- 
cles weights are updated to represent the new posterior. At each step in the process, a 
sum of the particles is taken to calculate a new state estimate. 

This study uses a standard particle filtering algorithm to estimate the state of the 
missile trajectory. The posterior distribution of the state is assumed to be Gaussian due to 
the additive Gaussian process noise. The steps to the algorithm are: 

1. Generate an initial set of N particles = 1,..., from the prior 

p(x 0 ) using Gaussian distribution about the estimated location of the mis- 
sile launch site with zero initial velocity and acceleration n each dimension. 

Each particle is given an equal normalized weight such that the i th particle 

weights is w' 0 = \/ N , where i=l,...,N. In this study N = 5000, or 10000. 
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2. For each particle, 

(a) Importance sampling step: Perform importance sampling on the parti- 
cles. A set of likelihood function l' k are generated from the difference be- 
tween the propagated particles and the observation. The corresponding 
weights are recursively updated based this likelihood function. 



(b) Selection step (resampling): Resampling of the new weights is accom- 
plished by comparing the particle weights to a threshold. Particles that fall 
below the threshold are discarded. Those that meet the criteria are resam- 
pled to keep the number of particles constant 

(c) The new state estimate is then calculated at each time by summing the 
particle values. Such that 



D. UNSCENTED PARTICLE FILTER (UPF) 

The unscented Kalman Filter is able to more accurately propagate the mean and 
covariance of the Gaussian approximation to the state distribution than the EKF. In com- 
parison to the EKF, the UKF tends to generate a more accurate estimate of the true co- 
variance of the state. Distributions generated by the UKF generally have a bigger support 
overlap with the true posterior distribution than the overlap achieved by the EKF esti- 
mates. This is in part related to the fact that the UKF calculates the posterior covariance 
accurately to the 3 rd order, whereas the EKF relies on a first order biased approximation. 
This makes the UKF a better candidate for more accurate proposal distribution generation 
within the particle filter framework. The UKF also has the ability to scale the approxima- 




( 11 ) 



N 




( 12 ) 



/= 1 
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tion errors in the higher tailed distributions. Because the sigma point set used in the UKF 
is deterministically designed to capture certain characteristics of the prior distribution, 
one can explicitly optimize the algorithm to work with distributions that have heavier 
tails than Gaussian distributions, i.e. Cauchy or Student-t distributions. This characteris- 
tic makes the UKF very attractive for the generation of proposal distributions [7], The 
algorithm of the Unscented Particle Filter is followings. 



1. Initialization: t = 0 

a. For i = 1, ..., N, draw the states (particle) x g from the prior /?(jc 0 )and 
set. 




P 0 (o 0 0 

0 Q 0 

0 0 7? 



2. For t = 1,2,... 

a. Importance sampling step 
• For i = 1, .., N: 

- Update the particle with the ULF: 

* Calculate sigma points: 

x±?=[*!T f'_ 

* Propagate particle into future (time update): 
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* Incorporate new observation (measurement update): 
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7=0 
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- Sample x, (,) □ q (x, (,) 1 xj',-, , , y u ) = N (*, (0 , ^ ) 

-Set4?0U--|.4 ) ) 4’(C,.4’) 

• For i =1, ... , N, evaluate the importance weights up to a normaliz- 
ing constant. 

b. Selection step 



• Multiply/Suppress particles (^o with high/low importance 
weights iv*' ) , respectively, to obtain N random particles (*ov ’-^a/*) 
c. MCMC step (optional) 
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• Apply a Markov transitional kernel with invariant distribution 
p(4! I y\,) to obtain {4U P Q<) • 

d. Output: The output is generated in the same manner as for the generic 
particle filter. 
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IV. SENSOR DATA FUSION 



Sensor data fusion is the process of combining outputs from sensors with informa- 
tion from other sensors, information processing blocks, database or knowledge bases, into 
one representational form. This technique is expected to achieve improved accuracy and 
more specific inferences than could be achieved by the use of a single sensor alone [14]. 
The applications of sensor data fusion include automated target recognition, guidance of 
autonomous vehicles, remote sensing, battlefield surveillance, automatic threat recogni- 
tion systems, monitoring of manufacturing processes, robotics and medical. The tech- 
niques employed in data fusion are drawn from diverse disciplines like digital signal 
processing, statistical estimation and control theory and classical numerical methods [14]. 
In principle, fusion of data from several sensors provides significant advantages over sin- 
gle source data. In addition to the statistical advantage gained by combining same-source 
data, the use of multiple types of sensors may increase the accuracy with which a quantity 
can be observed and characterized. Raw data from similar sensors may be directly com- 
bined if the sensors are measuring the same physical phenomenon. However feature/state 
vector or decision level fusion may be employed to combine data that come from non- 
commensurate sensors. 

While attempting to build a sensor data fusion system, the architecture that is cho- 
sen for fusion plays a significant role in deciding the algorithm for data fusion. The 
choice of architecture is made with a view to balancing computing resources, available 
communication bandwidth, desired accuracy and the capabilities of the sensors [14]. 

The three commonly used architecture are (i) centralized (ii) distrib- 
uted/decentralized and (iii) hybrid. The distributed/decentralized architectures have sev- 
eral inherent operational advantages that are dictated by the current trends towards modu- 
lar and autonomous systems [15]. For such architecture, information filtering which es- 
sentially tracks information about states is found to be most suitable. The information fil- 
ter is found to provide direct interpretation of node observations and contributions in 
terms of information. Prediction and estimation in terms of information reduce the com- 
putational load particularly for nonlinear systems in multi sensor systems. 
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The linear information filter (LIF) and extended information filter (EKF) are de- 
scried here. The information filter is a Kalman filter recast in terms of the information 
state vector and information matrix. The advantages of this formulation is that the esti- 
mation equations are computationally simpler than the measurement update equations of 
the Kalman filter, although prediction is more complex. Information filter is simpler to 
decouple and decentralize amongst a network of sensor nodes. A further advantage of 
the information filter is that it may be partitioned to provide a simple hierarchical estima- 
tion architecture based on the communication of the information terms from sensor nodes 
to a common fusion center and information state estimates from nodes to a common as- 
similation point. The proposed fusion architecture is also described in the end of this sec- 
tion. 



A. INTRODUCTION TO THE INFORMATION FILTER 

The diverse and sometimes conflicting information obtained from multiple sen- 
sors gives rise to the problem of how the information may be combined in a consistent 
and coherent manner. This is the data fusion problem. Multisensor fusion is the process 
by which information from multiple sensors is combined to yield a coherent description 
of the system under observation. A variety of information based data fusion algorithms 
have been employed in the recent work of Mutambara [12] and that of Mayika et al [13]. 

Most current sensor fusion algorithms consider systems described by linear dy- 
namics and observation models, whereas, most practical problems have nonlinear dynam- 
ics and sensor information nonlinearly dependent on the states that describe the environ- 
ment. Information-based estimation makes fully decentralized data fusion for nonlinear 
systems attainable. Consider a linear system with state-space representation, 

x(k) = F(k)x(k - 1) + B(k)u(k - 1) = w(k - 1) (l) 

Where %(k) is the state of interest at time k, F(k) is the state transition matrix form 
time (k-1) to k while, u(k) and B(k) are the input control vector and matrix, respectively. 
w(k) □ N(0,O(k)) is the associated process noise modeled as an uncorrelated, zero 
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mean and white sequence. The system is observed according to the linear discrete equa- 
tion 



z(k) = H (k)x(k) + v(k) (2) 

Where z(k) is the vector of observations made at time k. H(k) is the observation matrix or 
model and v(k) [ N(0,R(k )) is the associated observation noise modeled as an un- 
correlated white sequence. The estimate of the state x(j) at time i given information up to 
and including time j is given by 

x(i\j) = E[x(i)\z(\) i ...z(j)] 



p (i I j) = E (x(i)-x(i\ j)) (x(i)-x(i \ j)) T \ z(\),...z(j) 



B. LINEAR INFORMATION FILTER 

For a system described by equation (1) and being observed according to Equation 
(2), the Kalman filter provides a recursive estimate x(k | Ar ) for the state x(k) at time k 
given all information up to time k in terms of the predicted state x(k | k — 1) and the 

new observation z(k). The one-step-ahead prediction, x{k | k — 1), is the estimate of 
the state at a time k given only information up to time (k-1). The information filter is es- 
sentially a Kalman filter expressed in terms of measure of information about the parame- 
ters (state) of interest rather than direct state estimate and their associated covariances. 
This filter has also been called the inverse covariance form of the Kalman Filter. Assum- 
ing Gaussian noise and minimum mean squared error estimation gives the Fisher infor- 
mation matrix as, 

F{k) = p~\k\k) 

This information matrix or the inverse of the covariance matrix, is central to the filtering 
techniques. The general description of a linear system is given by: 
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X k + 1 - 0X k + w k 



( 3 ) 



Z k -HX k +v k ( 4 ) 

Where X k is the state vector at time k and <j) is the transition matrix. The input noise 
W* is assumed to be zero mean white Gaussian process with variance Q. The measure- 
ment noise P* j s assumed to be zero mean, white Gaussian with variance R. The in- 
formation state vector y(i|j) and the information matrix Y(ijj) are related to the state vec- 
tor and the covariance matrix P(i[j) by the equations: 

yd\j) = P~'a\j)xU\j) 

mj)=p-'a\j) 

Defining 

i(k) = H T (k)R~‘z(k) 

As the information state contribution from an observation z(k), and 

I{k) = H T {k)R-'H{k) 

As its associated information matrix and 

L(k\k-\) = P~' {k\k - 1 )F{k)P(k - 1 1 k - 1) 

As a propagation matrix (independent of the observation made), linear filter is 
written in terms of the information state/matrix: 



( 6 ) 



( 7 ) 



Prediction: 

y{k\k-\) = L(k\k-\)y(k-\\k-\) 
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( 11 ) 



Y(k\k-\) = F(k)Y~' (k - 1 1 k - \)F‘ (k) + Q(k ) 



Estimation: 



y(k\k) = y(k-\\k-\) + i(k) 



Y(k\k) = Y(k\k-\) + J(k) 



( 12 ) 

(13) 



The advantage i f this formulation is that the estimation equation (12) and (13) are compu- 
tationally simpler than the measurement update equations of the Kalman filter. 



C. EXTENDED INFORMATION FILTER 

The linear information filter can now be extended to a linearized estimation algo- 
rithm for nonlinear systems by using principles form both the derivations of the Informa- 
tion filter and the EKF. This generates a filter that predicts and estimates information 
about nonlinear state parameters given nonlinear observations and nonlinear system dy- 
namics. It is the Extended Information filter (EIF) and provides a solution to the nonlin- 
ear estimation problem. The EIF also has all the advantages of the Information filter and 
resolves some of the problems associated with the EKF. The estimation for nonlinear 
systems, in particular multisensory systems, is best carried out using information vari- 
ables rather than state variables. 

It is important to note that the EIF can’t be extrapolated from the Information fil- 
ter and EKF in an obvious manner. This is because in the nonlinear case, the function of 
the Information filter depends on this separation, which is possible in the linear observa- 
tion equation. In a real data fusion problem, the state of interest could evolve in a nonlin- 
ear manner, in which case simple linear models cannot adequately describe such a nonlin- 
ear system. Also, sensor observations may not be linear functions of the states. The gen- 
eral nonlinear system could be described by the following models. 
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X(k) = f(X(k-\),k) + w(k) 
Z(k) - h{X ( k ), k) + v(k) 



( 14 ) 



(15) 

Where X(k) is the state vector and w(k) is additive process noise vector at time step k. 
The nonlinear function f(.,k) is the nonlinear state transition function mapping the previ- 
ous state to the current state. Z is the observation vector at time k and h is the nonlinear 
observation model mapping the current state to the observations. The E1F algorithm is 
the followings: 

Prediction: 



y{k\k-\) = Y(k\k- \)f(k,ic(k - 



Y{k\k-\) = STf x (k)r\k-\\k- 



1 W/(k) + Q(k)V 



(16) 

(17) 



Estimation: 

y(k | k) = y(k - 1 1 k — 1 ) + i(k) 

Y(k\k) = Y(k\k-\) + I(k) 

The information state/matrix are given by 

I(k) = Vh T x (k)R~\k)Vh x (k) 
i{k ) = Vh T x ( k)R (*)[v(*) + Vh x (k)x(k |^-1)] 



(18) 

(19) 



( 20 ) 

( 21 ) 



where v(k) is the innovation sequence given by 

v(k) = z{k) - h(x (k | k - 1 )) (2 

The EIF described above has several attractive practical features: 

• The filter solves, in information space, the linear estimation problem for 

systems with both nonlinear dynamics and observations. In addition to having 
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all the attributes of the Information filter, it is a more practical and general filter. 

• The information estimation Equation (18) and (19) are computationally simpler 
Than the EKF estimation equations. This makes the partitioning of these 
equations for decentralized systems easy. 

• Although the EIF prediction Equation (16) and (17) are of the same apparent 
complexity as the EKF ones, they will be easier to distribute and fuse because of 
the orlhonormality properties of information space parameters. 

• Since the EIF is expressed in terms of information matrices and vectors, it is 
easily initialized compared to the EKF. Accurate initialization is important 
where linearized models are employed. 



Some of the drawbacks inherent in the EKF still affect the EIF. These include the 
nontrivial nature of Jacobian matrix derivation (and computation) and linearization insta- 
bility. 

D. DISTRIBUTED AND DECENTRALIZED DATA FUS101N SYSTEMS 

The nature of data fusion is that there are number of sensors physically distributed 
around an environment. In a centralized data fusion system, raw sensor information is 
communicated back to a central processor, where the information is combined to produce 
a single fused picture of the environment. In a distributed data fusion system, each sen- 
sor has it’s own local processor which can generally extract useful information from the 
raw sensor data prior to communication. This has the advantage that less information is 
normally communicated, the computational load on the central processor is reduced and 
the sensors themselves can be constructed in a reasonably modular manner. The degree 
to which local processing occurs at a sensor site varies substantially from simply valida- 
tion and data compression up to the full construction of tracks or interpretation of infor- 
mation locally. 
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While for many systems a centralized approach to data fusion is adequate, the in- 
creasing sophistication, functional requirements, complexity and size of data fusion sys- 
tems, coupled with the ever reducing cost of computing power argues more and more to- 
ward some form of distributed processing. The central issue in designing distributed data 
fusion systems is the development of appropriate algorithms which can operate at a num- 
ber of distributed sites in a consistent manner. 

This section begins with a general discussion of data fusion architectures and in- 
troduces the data fusion architecture used in this research. 

1. Data Fusion Architecture 

Distributed data fusion systems may take many forms. At the simplest level, sen- 
sors could communicate information directly to a central processor where it is combined. 
Little or no local processing of information need take place and the relative advantage of 
having many sources of information is sacrificed to having complete centralized control 
over the processing and interpretation of this information. As more processing occurs lo- 
cally, so computational and communication burden can be removed from the fusion cen- 
ter, but at the cost of reduced direct control of low-level sensor information. Increasing 
intelligence of local sensor nodes naturally results in a hierarchical structure for the fu- 
sion architecture. This has the advantage of imposing some order on the fusion process, 
but the disadvantage of placing a specific and often rigid structure on the fusion system. 
Other distributed architectures consider sensor nodes with significant local ability to gen- 
erate tracks and engage in fusion tasks. Fully decentralized architectures have no central 
processor and no common communication system. In such systems, nodes can operate in 
a fully autonomous manner, only coordinating through the autonomous communication 
information. 

In a hierarchical structure, the lowest level processing elements transmit informa- 
tion upwards, through successive levels, where the information is combined and refined, 
until at the top level some global view of the state of the system is made available. Such 
hierarchical structures are common in many organizations and have many well-known 
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advantages over fully centralized systems; particularly in reducing the load on a central- 
ized processor while maintaining strict control over sub-processor operations. 

The hierarchical approach to systems design has been employed in a number of 
data fusion systems and has resulted in a variety of useful algorithms for combining in- 
formation at different levels of a hierarchical structure. The Figures 1 1 and 12 show hi- 
erarchical fusion systems. A hierarchical approach to the design of data fusion systems 
also comes with a number of inherent disadvantages. The ultimate reliance on some cen- 
tral processor or controlling level within the hierarchical means that reliability and flexi- 
bility are often compromised. Failure of this central unit leads to failure of the whole 
system, changes in the system often mean changes in both the central unit and in all re- 
lated sub-units. Further, the burden placed on the central unit in terms of combining in- 
formation can often still be prohibitive and lead to an inability of the design methodology 
to be extended to incorporate an increasing number of sources of information. 









Sensors 




Figure 1 1 . Single Level Hierarchical Multiple Sensor Tracking System 
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Figure 12. Multiple Level Hierarchical Multiple Sensor Tracking System 



The move to more distributed, autonomous, organizations is clear in many infor- 
mation processing systems. This is most often motivated by two main considerations; the 
desire to make the system more modular and flexible, and a recognition that a centralized 
or hierarchical structure imposes unacceptable overheads on communication and central 
computation. The migration to distributed system organizations is most apparent in Arti- 
ficial Intelligence (AI) application areas, where distributed AI has become a research area 
in its own right. Many of the most interesting distributed processing organizations have 
originated in this area. 
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Figure 13. Blackboard Architecture in Data Fusion 

The above notable figure is the ‘'Blackboard” architecture, originally developed in 
the Hearsay speech understanding program, but now widely employed in many areas of 
AI and data fusion research. A Blackboard architecture consists of a number of inde- 
pendent autonomous “agents”. Each agent represents a source of expert knowledge or 
specified information facility or shared memory resource. This resource is called black- 
board. The blackboard is designed to closely replicate its physical analogue. Each agent 
is able to write information or local knowledge to this resource. Every agent in the sys- 
tem is able to read from this resource, in an unrestricted manner, any information which it 
considers useful in its current task. In principle, every agent can be made modular and 
new agents may be added to the system when needed without changing the underlying 
architecture or operation of the system as a whole. The flexibility of this approach to sys- 
tem organization has made the Blackboard architecture popular in a range of application 
domains. In data fusion, the Blackboard approach has been most widely used for knowl- 
edge-based data fusion systems in data interpretation and situation assessment. However, 
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the structured nature of tracking and identification problems does not lent itself to this 
anarchic organizational form. 

A decentralized data fusion system consists of a network of sensor nodes, each 
with its own processing facility, which together do not require any central fusion or cen- 
tral communication facility. In such a system, fusion occurs locally at each node on the 
basis of local observations and the information communicated from neighboring nodes. 
At no point is there a common place where fusion of global decisions are made. 

A decentralized data fusion system is characterized by three constraints: 

• There is no simple single central fusion center; no one should be central to the 
successful operation of the network. 

• There is no common communication facility; nodes cannot broadcast results and 
communication must be kept on a strictly node-to-node basis. 

• Sensor nodes do not have any global knowledge of sensor network topology; 
nodes should only know about connections in their own neighborhood. 



Figure 14 show three possible realizations of a decentralized data fusion system. 
The key point is that all these systems have no central fusion center. 

The constraints imposed provide a number of important characteristics for decen- 
tralized data fusion systems: 

• Eliminating the central fusion center and any common communication facility 
ensures that the system is scalable as there are no limits imposed by centralized 
computational bottlenecks or lack of communication bandwidth. 

• Ensuring that no node is central and that no global knowledge of the network 
topology is required for fusion means that the system can be made survivable to 
the on-line loss (or addition) of sensing nodes and to dynamic changes in the 
network structure. 
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• As all fusion processes must take place locally at each sensor site and no global 
knowledge of the network is required a priori, nodes can be constructed and 
programmed in a modular fashion. 

These characteristics give decentralized systems a major advantage over more 
traditional sensing architectures, particularly in defense applications. 




Figure 14. Blackboard Architecture in Data Fusion 
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Figure 15. Blackboard Architecture in Data Fusion 




Figure 16. Blackboard Architecture in Data Fusion 
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2. Proposed Fusion Architecture 

The proposed hierarchical fusion architecture is a simple hierarchical estimation 
architecture based first on the communication of the information terms /(□) and /(C) 

from sensor nodes to a common fusion center. In this system, information-state contribu- 
tions are calculated at each sensor node and transmitted to a central fusion center where a 
common estimate is obtained by simple summation. All state predictions are undertaken 
at the central processor. Each sensor incorporates a full state model and takes observa- 
tions according to the Equation (2). They all calculate an information state contribution 
from their observations in terms of /, (&) and /, (£) . These are then communicated to the 

fusion center and are incorporated into global estimate through Equations (18) and (19). 
The information-state prediction is generated centrally using Equations (16) and (17) and 
the state estimate itself may be found at any stage from Jc(/ 1 j) = Y~' (/| y) j>(/| /) . To 

avoid communicating predictions to nodes, any validation or data association should take 
place at the fusion center. Figure 17 shows the hierarchical architecture. 



z' '■ 




Sensor 



Figure 17. A Proposed simple hierarchical fusion architecture 
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V. SIMULATION RESULTS 



In this section, we will describe the results of simulation studies on the perform- 
ance of nonlinear filter and sensor fusion using the proposed fusion architecture. The 
IMPULSE model and MATLAB are the tools used to implement the simulation. Specifi- 
cally, the results from the IMPULSE program are used by the MATLAB code from [16] 
in order to simulate the tracking problem. The number of trajectories used in the simula- 
tion is totally 12 including upper and lower stage for each of 6 ballistic missile trajecto- 
ries. A detailed description of the requirements (inputs) and capabilities (outputs) of the 
IMPULSE program are outlined in [16]. An unclassified model of the ballistic missile is 
used. Figure 23 shows the trajectories. 




Figure 1 8. Six near-simultaneous launches of ballistic missiles (viewed looking north). 
The delay ranges from four to ten second intervals. Missile staging occurs at 65 to 85 
seconds, simulation time. The coordinate system here is local vertical (north-east-up) 
with respect to the sensor. 
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Figure 19. Same profile as in Figure 23: Six, near simultaneous, ballistic missile launches 
(looking East-to-West) and 7 RF sensors 



Table 4 lists the positions of the sensors in the vicinity of the launched ballistic 
missile. The radars are placed at sea level. 



Sensors 


Position 


Latitude 


Longitude 


Altitude 


RF #1 


131°46'44 "E 


40°50'45 "N 


Sea level 


RF #2 


132°46'44"E 


41 o 40'45"N 


Sea level 


RF #3 


1 3 TOO '00 "E 


40°00'00"N 


Sea level 



45 



RF #4 


130°00'00"E 


35°00'00"N 


Sea level 


RF #5 


135°00'00" E 


40°00'00"N 


Sea level 


RF #6 


140°00'00"E 


45°00'00"N 


Sea level 


RF #7 


145°00'00"E 


50°00'00"N 


Sea level 



Table 4. Sensor positions 



A. MODELS OF TARGET MOTION AND RADAR MEASUREMENTS 

A ballistic missile tracking algorithm is EIF based on the extended Kalman form. 
The system equations are the standard tracking equations 

x M =F k x k +G k +co k 
h = K x u + v k 

where the missile state vector is given as 
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The term F h is the linear state transition matrix: 
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The gravity matrix, G k , which accounts for the acceleration in z due to gravity 



with £ =9.8, m ■ s ' 1 , is 






~g* 
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0 

0 

0 

A^_ 

2 

A 

0 



and co k is the plant noise with covariance O k 
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where A is the sampling interval and q 2 is a scaling factor used to account for un- 
modeled target maneuver accelerations, and v k is the measurement noise with covariance 



R, = 



a R 0 



0 

0 



0 

0 



Although the missile dynamics in this system are linear, the measurement process 
is nonlinear. The sensor platform observes the missile positions through measurements 
in range, azimuth and elevation relative to the sensor. 



Let 
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K = 



Range 

Azimuth 

Elevation 



where 



range 



= Jx 2 +y 2 +z 2 = + JV4 + -7 



azimuth = tan 1 


V 


= tan -1 


*4 




x _ 







with 



elevation = tan 




r' = Jx 2 +y 2 = Jxf +yt 

These measurement equations are nonlinear and therefore must be converted us- 
ing a series expansion of the measurement equation h k . Applying the definition of the 

H r k matrix, the gradient of h k is determined as follows: 



dr(x) 


8r(x) 


dr(x) 


dr( x) 


dr(x) 
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dr(x) 
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5x 9 
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dx 9 



which simplifies to 
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Finally, the approximate — linearized — system of equations for the state predict and mea- 
surement are the followings. 



**♦1 = F k*k +G k+“>k 



z k = #*** +»'* 



b.comparison of performance of nonlinear filters 

In this section, we describe the experimental results of comparing the perform- 
ance of three nonlinear Filters, EKF, UKS, and UPF , since the PF did not converge. We 
used the distance error, which is the square root of sum for each coordinate distance er- 
ror, as a measure to evaluate the performance. In this study, two values for the measure- 
ment noise are chosen to have the following standard deviations: 

^ range ^azimuth ® elevation ^ 



Grange = 1 0 meters, cr^'b = T 



10 

u elevation 



Figures 21-24 show the average distance errors for two of the trajectories in the 
experiments over 30 Monte Carlo runs. Each sensor performs the estimation for each tra- 
jectory and then the results are averaged for the time period. The red line shows the dis- 
tance error calculated from each coordinates after EKF estimation is done. The green 
line shows the results for the UKF and the blue line is for the UPF. The diamond marker 
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is designating the distance error at the last time. As the graphs show, when standard de- 
viations for measurement noise are all set to 0, the EKF is the best. 



51 



Figure 20. 



Figure 2 1 . 



Average Disance error of 7 sensors for U1 trajectory 




Average distance error of seven sensors for 1 st upper stage trajectory 

( ^ range ^ azimuth ^ elevation ^ ) 



Average Distance error of 7 sensors for LI trajectory 




Average distance error of seven sensors for 1 st lower stage trajectory 

( ^ range ^ azimuth ^ elevation ^ ) 
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Figure 22. 



Figure 23. 



Average Disance error of 7 sensors for U2 trajectory 




Average distance error of seven sensors for 2 nd upper stage trajectory 

( ^ range ^"azimuth ^ elevation ^ ) 



Average Distance error of 7 sensors for L2 trajectory 




Average distance error of seven sensors for 2 nd lower stage trajectory 

( ® range ^ azimuth ^ elevation ^ ) 
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Table 5 shows average distance error and final distance error at last time for all 
trajectories in number 



Trajectory 


EKF 


UKF 


UPF 


Average 


Last 


Average 


Last 


Average 


Last 


U1 


0.0089 


0.0238 


0.0258 


0.0237 


0.0223 


0.0399 


LI 


0.0034 


0.0001 


0.0308 


0.0264 


0.0222 


0.0203 


U2 


0.063 


0.0381 


0.0281 


0.0506 


0.0222 


0.0334 


L2 


0.0026 


0.0007 


0.037 


0.028 


0.0247 


0.0226 


U3 


0.0076 


0.0231 


0.0259 


0.0226 


0.0214 


0.035 


L3 


0.0036 


0.0008 


0.0323 


0.0283 


0.0229 


0.0152 


U4 


0.0166 


0.2043 


0.0315 


0.0988 


0.0292 


0.0938 


L4 


0.0013 


0.0007 


0.029 


0.0212 


0.0201 


0.0142 


U5 


0.0079 


0.0559 


0.037 


0.0316 


0.0273 


0.0259 


L5 


0.0047 


0.0006 


0.0374 


0.033 


0.026 


0.0175 


U6 


0.0067 


0.0078 


0.0282 


0.0595 


0.0225 


0.041 


L6 


0.0037 


0.0006 


0.0385 


0.0325 


0.0261 


0.0253 



Table 5. Average distance error and distance error at last time of three nonlinear filters 
for 1 2 trajectories ( cr rmge = <J mimu , h = a eMon = 0 ) 

However, the performance of the EKF is degraded very badly for the second case. 
Figures 24-27 shows the average distance error for four of the trajectories. The second 
figure compares the distance error of only UPF and UPF. They work well without the 
performance degradation for the second case, as the results show. Even though the UPF 
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has the best performance of three filters, because of its time complexity, we can know 
that the UKF is the best estimation algorithm for ballistic missile tracking through . 



Average Distance error of 7 sensors for LI trajectory 




Figure 24. Average distance error of seven sensors for 1 st lower trajectory 

(°W =10 meters, a a:wwlh = 1 V c/mM(0 „ = 1°) 
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Average Distance error of 7 sensors for LI trajectory UKF v.s UPF 




Figure 25. Distance error for the 1 sl lower stage trajectory (UKF vs. UPF) 

( Grange = 1 0 meters, a a:mmlh = 1°, a clevollon = 1° ) 



Average Disance error of7 sensors for U1 trajectory 




Figure 26. 



Average distance error of seven sensors for 1 upper trajectory 

( ^ range = 1 = lV„„ = 1°) 
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Average Distance error of 7 sensors for U1 trajectory. UKF v.s UPF 




Figure 27. Distance error for the 1 st upper stage trajectory (UKF vs. UPF) 

( ^ range = 10 meters, a a:imulh = 1 = 1°) 



Trajectory 


EKF 


UKF 


UPF 


Average 


Last 


Average 


Last 


Average 


Last 


U1 


16.3177 


89.0724 


0.0244 


0.0368 


0.0149 


0.0364 


LI 


10.2068 


0.3315 


0.0311 


0.0259 


0.0147 


0.0089 


U2 


18.9756 


160.6896 


0.0281 


0.0765 


0.0182 


0.0713 


L2 


25.2276 


0.0115 


0.1107 


0.028 


0.1007 


0.026 


U3 


16.111 


48.2313 


0.0258 


0.0301 


0.0167 


0.0447 


L3 


9.0777 


3.0363 


0.0324 


0.0298 


0.0168 


0.0116 


U4 


28.1213 


179.696 


0.0339 


0.1223 


0.0197 


0.0163 
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L4 


2.6186 


0.0115 


0.0287 


0.0214 


0.0141 


0.0055 


U5 


12.4995 


36.0934 


0.0302 


0.0542 


0.0206 


0.0261 


L5 


15.1794 


1.1 139 


0.0368 


0.0332 


0.0196 


0.0251 


U6 


16.8421 


38.389 


0.0284 


0.0353 


0.0173 


0.0141 


L6 


18.1476 


7.1747 


0.0385 


0.0325 


0.0206 


0.0095 



Table 6. Average distance error during the total time interval and average distance error 
at last time of three nonlinear filters for 12 trajectories 

( G range ~ ^ OtH, & — 1 , G e i erannn = 1 ) 



C. MULTI-SENSOR FUSION 



We compare the performance of the sensor fusion method using the EIF algo- 
rithm. We used the distance error, which is the square root of the sum for each coordinate 
distance error, as a measure to evaluate the performance. In this study, the measurement 
noise is chosen to have the following standard deviation: 



■ cr 



range 



= 1 0 meters 



^ azimuth ^ 



^ 'elevation ^ 



In order to evaluate the performance of the sensor fusion method, the Track Score 
Function (TSF) proposed in the previous master’s theses is considered. The likelihood 
ratio for a combination of data, taking the priori probability data into account, is given as 
[19]: 



a (H v H,) = 



P{D!H y )Po{H x ) „ 
P(D/H 0 )Po(H 0 ) 



Ft 

Pf 
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Where //, is the true target hypothesis, H 0 is the false target hypothesis, P T is the prob- 
ability of true target, P r is probability of false target, P(D\H X ) is the probability density 
function assuming that the //, is correct, Po{H x ) is priori probability of //, . 

The performance of the tracking system should be independent of the behavior of 
the target. The system is evaluated for its ability to respond instantly to any change of the 
trajectory of the target either in velocity or direction. 

Taking K scans of data, where the measurement error for one is not related to the 
error from the previous scan, the likelihood ratio (a) can be given as the product of the 
likelihood ratio due to the kinematics (a a ) and the likelihood ratio due to the signal (a,) 
is given: 

K 

a (^) = a 0 []a < a 5 

K= 1 

Where 



PAH,) 

PAH o) 



In this experiment, we just considered a signal-related likelihood ratio. Based on 
likelihood ratio of the signal-related data, the track score is given by [1 8, Chapter 6]: 



A PjDel/HMyJDel-H,) 

■ P{Det / H 0 )p{y s l Det,H 0 ) 



(I) 



Taking both P D and P FA into account, (1) can be written as: 



A _ PpPiys 1 Det ^H x ) 

P F ,p(yJ DetM 0 ) 

In the case of radar, when the SNR is available, the signal likelihood ratio is given 
by (1). The probability density function under //, us given by 
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p(y/H i) 



1+v/ (1 + 2/0) f —v 
— ^-exp — = — 

(1 + 0/2) 2 P Ll + ©/2j 



Where y is SNR data and © is average of the SNR. 

Probability density function under H 0 us given by 

P(y/H 0 ) = exp(-.y) 



Figures 25-47 show the distance error in the experiments that completed 30 
Monte Carlo runs. The yellow line shows the distance error calculated from the sensor 
with the smallest error for three coordinates, where each sensor tracks the same target re- 
spectively. The blue line and green line are related to the TSF. The one results from a 
sensor with the biggest TSF at each time interval. The other one shows the results that 
fuse the top four sensors according to the TSF. Finally, the red line fuses all seven sen- 
sors. Figure 26 shows only the results of fusing all sensors and four sensors. Two fig- 
ures of results are generated for each trajectory. The first figure shows the distance error 
graph of three fusion method and that of sensor with the smallest error at each time inter- 
val. The second one is to enlarge the result of only all sensor and four sensor fusion 
method to show the distance error. The diamond marker is designating the distance error 
at the final time. 
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Distance error for U1 tracjectory 




Figure 28. Distance error for the 1 st upper stage trajectory (Three fusion methods and 
Sensor with smallest error) 



Distance error for U1 tracjectory 




Figure 29. Distance error for the 1 st upper stage trajectory (All sensors vs. Four sensors) 
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Distance error for LI tracjectory 




Figure 30. Distance error for the 1 st lower stage trajectory (Three fusion methods and 
Sensor with smallest error) 

Distance error for LI tracjectory 




Figure 3 1 . Distance error for the 1 st lower stage trajectory (All sensors vs. Four sensors) 
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Distance error for U2 tracjectory 




Figure 32. Distance error for the 2 nd upper stage trajectory (Three fusion methods and 
Sensor with smallest error) 



Distance error for U2 tracjectory 




Figure 33. Distance error for the 2 nd upper stage trajectory (All sensors vs. Four sensors) 
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Distance error for l_2 tracjectory 




Figure 34. Distance error for the 2 nd lower stage trajectory (Three fusion methods and 
Sensor with smallest error) 

Distance error for 1_2 tracjectory 




Figure 35. Distance error for the 2 nd lower stage trajectory (All sensors vs. Four sensors) 
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Distance error for U3 Iracjectory 




Figure 36. Distance error for the 3 rd upper stage trajectory (Three fusion methods and 
Sensor with smallest error) 



Distance error for U3 tracjectory 




Figure 37. Distance error for the 3 rd upper stage trajectory (All sensors vs. Four sensors) 



65 



Distance error for L3 tracjectory 




Figure 38. Distance error for the 3 rd lower stage trajectory (Three fusion methods and 
Sensor with smallest error) 



Distance error for LB tracjectory 




Figure 39. Distance error for the 3 rd lower stage trajectory (All sensors vs. Four sensors) 
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Distance error for U4 tracjectory 




Figure 40. Distance error for the 4 th upper stage trajectory (Three fusion methods and 
Sensor with smallest error) 



Distance error for U4 tracjectory 




Figure 41 . Distance error for the 4 th upper stage trajectory (All sensors vs. Four sensors) 
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Distance error for L4 tracjectory 




Figure 42. Distance error for the 4 th lower stage trajectory (Three fusion methods and 
Sensor with smallest error) 



Distance error for L4 tracjectory 




Figure 43. Distance error for the 4 th lower stage trajectory (All sensors vs. Four sensors) 
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Distance error for U5 tracjectory 




Figure 44. Distance error for the 5 th upper stage trajectory (Three fusion methods and 
Sensor with smallest error) 



0.25 



0.15 - v 






005 



Distance error for U5 tracjectory 




Figure 45. Distance error for the 5 th upper stage trajectory (All sensors vs. Four sensors) 
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Distance error for L5 tracjectory 




Figure 46. Distance error for the 5 th lower stage trajectory (Three fusion methods and 
Sensor with smallest error) 



Distance error for L5 tracjectory 




Figure 47. Distance error for the 5 th lower stage trajectory (All sensors vs. Four sensors) 
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Distance error for U6 tracjectory 




Figure 48. Distance error for the 6 th upper stage trajectory (Three fusion methods and 
Sensor with smallest error) 



Distance error for U6 tracjectory 




Figure 49. Distance error for the 6 lh upper stage trajectory (All sensors vs. Four sensors) 



71 



Distance error for U6 tracjectory 




Figure 50. Distance error for the 6 th lower stage trajectory (Three fusion methods and 
Sensor with smallest error) 



Distance error for U6 tracjectory 




Distance error for the 6 th lower stage trajectory (All sensors vs. Four sensors) 
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Figure 5 1 . 



Tables 7 and 8 shows average distance error during total time period for three sen- 
sor fusion methods and sensor with the smallest error at each time interval. The result of 
fusion of all sensors or four sensors is better than the result of sensor with the best per- 
formance. Also, we can see that the performance of all sensor fusion is only slightly bet- 
ter than that of four sensor fusion. In conclusion, we can get good performance through 
sensor fusion based on the Extended Information Filter and can overcome the drawbacks 
of the Extended Kalman Filter in terms of performance and execution time. 



Trajectory 


All sensors 


Four sensors 


TSF (Track 
Score Function) 


Sensor with the 
smallest error 


U1 


0.0225 


0.0238 


3.2469 


1.4361 


U2 


0.0385 


0.0412 


3.5311 


1.5719 


U3 


0.0406 


0.0437 


3.1981 


1.5035 


U4 


0.0294 


0.0311 


3.9471 


1.532 


U5 


0.0887 


0.0913 


5.1942 


1.8709 


U6 


0.0177 


0.0184 


3.4899 


1.5368 


LI 


0.0444 


0.0468 


3.8185 


1.6734 


L2 


0.0644 


0.0681 


4.327 


1 .8802 


L3 


0.0986 


0.0996 


4.9319 


1.6297 


L4 


0.0208 


0.0218 


3.5549 


1.6639 


L5 


0.0359 


0.0425 


3.3999 


1.4595 


L6 


0.0625 


0.0598 


4.4784 


1.8106 



Table 7. Average distance error of Monte-Carlo runs for 12 trajectories (Km) 
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Trajectory 


All sensors 


Four sensors 


TSF (Track 
Score Function) 


Sensor with the 
smallest error 


U1 


0.0165 


0.0302 


6.1695 


3.0208 


U2 


0.0434 


0.0667 


6.7872 


2.8399 


U3 


0.0407 


0.0397 


6.7577 


3.3037 


U4 


0.0297 


0.03087 


6.5375 


3.1721 


U5 


0.0915 


0.1157 


10.9499 


3.5601 


U6 


0.0166 


0.0232 


6.2319 


2.8892 


LI 


0.0273 


0.0294 


6.3306 


2.8978 


L2 


0.033 


0.0439 


7.7776 


3.4658 


L3 


0.0581 


0.1807 


10.1092 


3.4979 


L4 


0.0104 


0.0237 


5.4433 


3.3111 


L5 


0.0393 


0.0281 


5.9995 


2.9402 


L6 


0.056 


0.1 149 


8.249 


3.4377 



Table 8. Average distance error at the last time of Monte-Carlo runs for 12 trajectories 

(Km) 
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VI. CONCLUSION 



This study has evaluated the performance of nonlinear filters, the Extended Kal- 
man Filter, the Unscented Kalman Filter, the Particle Filter, and the Unscented Particle 
Filter and also developed a sensor fusion method to track multi-ballistic targets during 
boost phase using multi-sensors. This research makes use of the IMPULSE© simulation 
tool for generating ballistic missile flight profiles which serve as the test data for the 
tracking. We use the distance error as a performance measure. Our experimental results 
of comparing the performance of three nonlinear filters shows that the UPF has the best 
performance, but that the UKF is the best choice in tenns of time, complexity, and per- 
formance, while the standard PF does not converge. 

The fusion of data from several sensors provides significant advantages over sin- 
gle source data. In addition to the statistical advantage gained by combining same-source 
data, the use of multiple types of sensors may increase the accuracy with which a quantity 
can be observed and characterized. In this study, seven active ground-based radar sen- 
sors are used to track the ballistic missile during boost phase. We consider the Informa- 
tion Filter which is an efficient form to fuse the information from multiple sensors with- 
out information loss. To evaluate the performance of the Extended Information Filter, we 
compare with the track score function proposed in Patsikas’s previous work in which the 
calculation of the track scoring function is to identify the sensor with the best track file. 
This study tells us that the result of multiple sensor fusion outperforms the performance 
of one sensor with the best performance and the information form of filter has an advan- 
tage in constructing a more efficient fusion architecture which is hierarchical and easily 
lends itself to decentralized processing. 
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VII. FUTURE WORK 



In future work, the fusion architecture needs to be extended to a more efficient, 
fully decentralized form, so we may evaluate its performance and efficiency with those of 
the hierarchical fusion architecture. In addition, we should apply the Information Filter 
based on the Unscented Kalman Filter to track the Boost-phase ballistic missile, instead 
of the Extended Information Filter. 
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